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ABSTRACT 

We present three-dimensional numerical simulations of the rise and fragmentation of 
twisted, initially horizontal magnetic flux tubes which evolve into emerging fi-loops. The 
flux tubes rise buoyantly through an adiabatically stratified plasma that represents the 
solar convection zone. The MHD equations are solved in the anelastic approximation, 
and the results are compared with studies of flux tube fragmentation in two dimensions. 
We find that if the initial amount of field line twist is below a critical value, the degree of 
fragmentation at the apex of a rising il-loop depends on its three-dimensional geometry: 
the greater the apex curvature of a given O-loop, the lesser the degree of fragmentation 
of the loop as it approaches the photosphere. Thus, the amount of initial twist necessary 
for the loop to retain its cohesion can be reduced substantially from the two-dimensional 
limit. The simulations also suggest that as a fragmented flux tube emerges through a 
relatively quiet portion of the solar disk, extended crescent-shaped magnetic features 
of opposite polarity should form and steadily recede from one another. These features 
eventually coalesce after the fragmented portion of the Jl-loop emerges through the 
photosphere. 

Subject headings: methods: numerical — MHD, Sun: interior. Sun: magnetic fields 

1. Introduction 

The largest concentrations of magnetic flux on the Sun occur in active regions. Great progress 
has been made over the past decade in understanding the connections between the magnetic field 
in active regions, observed at the surface of the Sun, to the magnetic field deep in the solar interior. 
Active regions have a bipolar structure, suggesting that they are the tops of magnetic flux loops 
which have risen from deep in the solar interior. On average, active regions are oriented in the 
E-W direction (Hale's Polarity Law) suggesting that the underlying field geometry is toroidal. 
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The persistence of Hale's law for periods of several years during a given solar cycle suggests that 
magnetic flux must be stored in a relatively stable region of the solar interior. Several stability 
arguments (Spruit & van Ballegooijen 1982; van Ballegooijen 1982; Ferriz-Mas & Schiissler 1993, 
1995) show that the only place where such fields can be confined stably for periods of several years 
is below the solar convection zone. On the other hand, if magnetic fields are placed any significant 
distance below the top of the radiative zone, they are so stable they could not emerge on the time 
scale of a solar cycle. We are thus led to the conclusion that the most likely origin of active region 
magnetic fields is from a toroidally oriented field layer residing in the "convective overshoot region" , 
a thin, slightly convectively stable layer just beneath the convection zone. This layer also seems to 
coincide with the "tachocline" (Kosovichev 1996; Corbard et al. 1999), where the solar rotation rate 
transitions from solid body behavior in the radiative zone, to the observed latitudinally dependent 
rotation rate we see at the Sun's surface. This suggests that not only are solar magnetic fields 
stored in the convective overshoot layer, the overshoot layer is also the most likely site for the solar 
cycle dynamo (Oilman, Morrow, & DeLuca 1989; DeLuca &: Oilman 1991; Parker 1993; MacGregor 
&; Charbonneau 1997a, b; Durney 1997; Dikpati & Charbonneau 1999). 

Over the past decade, most efforts to study the emergence of active region magnetic fields 
have employed the "thin flux tube" approximation. This model assumes that magnetic flux tubes 
behave as distinct tube-like entities, surrounded by field-free plasma. The approximation further 
assumes that the tube diameter is small compared to all other length scales in the problem, and 
that pressure balance exists across the tube at all times. After adopting these assumptions, it is 
straightforward to derive an equation of motion for the dynamics of the tube from the momentum 
equation in MHD. Thin flux tube models of emerging active regions have proven very successful in 
explaining many properties of active regions in terms of flux tube dynamics in the solar interior. For 
example, they have successfully explained the variation of active region tilt with respect to the E-W 
direction as a function of solar latitude (D'Silva & Howard 1993; D'Silva k, Choudhuri 1993; Fisher, 
Fan k, Howard 1995), the asymmetric orientation of the magnetic field after emergence (van Driel- 
Gesztelyi & Petrovay 1990; Moreno-Insertis, Schiissler, k, Caligari 1994; Cauzzi, Moreno-Insertis, 
& van Driel-Oesztelyi 1996), and the observed scatter in tilts as a function of active region size 
(Longcope & Fisher 1996). 

In spite of these successes, recent two-dimensional MHD simulations of flux tube emergence 
have shown results which seem to invalidate many assumptions that are adopted in the thin flux 
tube approximation. Schiissler (1979) and Longcope, Fisher, k. Arendt (1996) find that an initially 
buoyant, untwisted flux tube will fragment into two counter-rotating tube elements which then 
separate from one another, essentially destroying the tube's initial identity. Moreno-Insertis k 
Emonet (1996) and Fan, Zweibcl, k Lantz (1998) have demonstrated via two-dimensional MHD 
simulations that in order to prevent a flux tube from fragmenting, enough twist must be introduced 
into the tube to provide a cohesive force to balance the hydrodynamic forces acting to rip it 
apart. That critical twist is defined, roughly, by that necessary to make the Alfven speed from the 
azimuthal component of the field at least as great as the relative velocity between the tube and 
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the field free plasma surrounding it (Linton, Longcope, & Fisher 1996; Emonet &; Moreno-Insertis 
1998; Fan, Zweibel, & Lantz 1998). 

But when global levels of twist in active regions (Pevtsov, Canfield, k. Metcalf 1995; Longcope, 
Fisher, & Pevtsov 1998) are determined from vector magnetograms, the amplitude of the observed 
twist is typically far smaller than this critical value (Longcope et al. 1999). When plotted as 
a function of active region latitude, the twists exhibit large scatter, but superimposed on this 
apparently random behavior there is a slight, but clearly discernible trend for active regions in the 
northern hemisphere to be negatively twisted, while those in the south are positively twisted. 

Longcope, Fisher, & Pevtsov (1998) have developed a theoretical model which not only explains 
the latitudinal variation of twist, but also can account for the large fluctuations in twist from active 
region to active region. In this model, an initially untwisted flux tube rises through the convection 
zone in accordance with the thin flux tube approximation. Coriolis forces acting on convcctive 
eddies produce a non-zero average kinetic hclicity, which is proportional to latitude. The kinetic 
helicity acts to "writhe" the flux tube, which is then twisted in the opposite direction to preserve 
its magnetic helicity. Longcope, Fisher, k, Pevtsov (1998) showed that this model can explain the 
observed data. Yet the model assumes from the beginning that the thin flux tube approximation 
can be used, even for an initially untwisted tube, while the two-dimensional MHD simulations 
suggest that this is invalid. Is there some way out of this quandary, which we dub "Longcope's 
Paradox" ? 

In this paper, we describe MHD simulations of flux tube fragmentation in three dimensions. 
The result of these simulations is that the critical degree of twist necessary to prevent fragmentation 
is reduced dramatically by the presence of flux tube curvature, as will be present in an emerging $7- 
loop. We find that for a fixed amount of twist, the degree of fragmentation is a function of the tube's 
curvature, and transitions asymptotically to the two-dimensional limit as the curvature approaches 
zero. Even for fiux tubes with little to no initial twist, the fragmenting magnetic morphology at 
the apex of an emerging Q-loop is considerably less dispersed than the two-dimensional simulations 
would indicate — a finding consistent with the simulation results of Dorch &: Nordlund (1998). It 
is not clear at this time whether our results resolve Longcope's paradox or not, but they certainly 
ameliorate the problem a great deal. 

The remainder of this paper is organized as follows: In Section 2 we briefly discuss the formalism 
of the anelastic approximation employed in our models, together with the numerical methods used 
to solve the system of equations. We also describe the range of initial configurations that we use 
to explore the relationship between flux tube fragmentation, tube geometry, and the initial twist of 
the fleld lines. At the beginning of Section 3, we define what is meant by a flux tube in the context 
of our three-dimensional MHD simulations, and further define a quantitative measure of the degree 
of fragmentation of such a tube. We then present the results of our numerical simulations and 
discuss the implications of the models. Finally, in Section 4 we summarize our conclusions. 
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2. Method 

We solve the three-dimensional MHD equations in the anelastic approximation (Ogura & 
Phillips 1962; Gough 1969) using a portable, modularized version of the code of Fan et al. (1999). 
The anelastic equations result from a scaled-variable expansion of the equations of compressible 
MHD (for details see Lantz & Fan 1999 and references therein), and describe variations of a plasma 
about a stratified, isentropic reference state. This approach is valid for low acoustic Mach number 
plasma below the photosphere (M = v/cs <C 1), where the Alfvcn speed fa is much less than the 
local sound speed Cg (Glatzmeier 1984). The primary computational advantage of this technique 
is that fast-moving acoustic waves are effectively filtered out of the simulations. This allows for 
much larger timesteps than would be possible in fully compressible MHD, and unlike a Bousinnesq 
treatment, a non-trivial background stratification can be included in the models. This method is 
well suited to our investigation of the evolution and fragmentation of magnetic flux tubes, since our 
region of interest lies well within the limits of validity of this approximation, and since we require 
many simulations to fully explore the relevant parameter space. It is important to note, however, 
that where M is not ^ 1, a fully compressible treatment (such as that of Bercik, Stein, &; Nordlund 
1999) is required. The equations of anelastic MHD are as follows: 

V • (pov) = (1) 
+ V • Vv j = -Vpi + pig + — (V X B) X B + V • n (2) 

poTo + V • V (so + si)^ = V • (KpoToVsi) + |- |V x B|2 + (n ■ V) • v (3) 

V • B = (4) 

5B 

_ = V X (v X B) V X (jyV X B) (5) 

Pl = Pl-Il (6) 
Po Po Tq 

Cp To 1 Po' 

Here, pi, pi, Ti, si, v, and B refer to the density, gas pressure, temperature, entropy, velocity, and 
magnetic field perturbations, while po, po, Tq, and sq denote the corresponding values of the zeroth 
order reference state, described in detail below, g = — g z is the acceleration due to gravity, and is 
assumed to be uniform in our calculations. The quantity Cp represents the specific heat at constant 
pressure. The viscous stress tensor 11 is given by 

and /X, 77, K, represent the coefficients of viscosity, magnetic diffusion, and thermal diffusion re- 
spectively. 
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A detailed description of the numerical methodology we use in the solution of the anelastic 
MHD equations can be found in Appendix A of Fan et al. (1999). Briefly, the non-dimensional 
form of the equations are solved in a rectangular domain assuming periodic boundary conditions 
in the horizontal directions, and non-penetrating, stress-free conditions at the upper and lower 
boundaries. The magnetic field B and the momentum density are both divergence-free, and 
thus each can be expressed in terms of two scalar potentials: 



These potentials, along with the other dependent variables of the problem, are spectrally decom- 
posed in the horizontal Cartesian directions. The Fourier variables are discretized with respect to 
the vertical direction, and the vertical derivatives are approximated by fourth-order, centered dif- 
ferences. A semi-implicit method is then used to time-advance the five discretized scalar equations 
for the Fourier variables. Using operator splitting, the second-order Adams-Bashforth scheme is 
applied to the advcction terms, and the second-order Crank-Nicholson scheme is applied to the 
diffusion terms (sec Press et al. 1986 for a general discussion of these methods). 

To investigate the dynamics of flux tube fragmentation, and how this process depends on the 
initial state and eventual geometry of the tube, we carried out a total of 16 simulations. In each case, 
an ideal gas of 7 = 5/3 is assumed; the reference state is taken to be an adiabatically stratified 
polytrope of index m = 1.5 (related to 7 by m = 1/(7 — 1)); and /x, rj, and poK are assumed 
constant throughout the simulation domain. The diffusive parameters enter into the calculation 
via the Reynolds number (i?c = Wi^]^]/^)) the magnetic Reynolds number {Rm = [-2] /??) , and 
the Prandtl number (Pj. = h/{Kpq)). The density and temperature scales {[p] and [T]) are defined 
as the density and temperature of the reference state at the bottom of the simulation domain (pr 
and Tr), and the length scale [z] is defined as the pressure scale height of the reference state at that 
same location (Hr = r^fTj-jg, where r* = R/li, and R and 11 are the ideal gas constant and mean 
molecular weight respectively). The velocity scaling [v\ is given as the characteristic Alfven speed 
along the axis of the initial magnetic flux tube. For each simulation, both R^ and i?ni are set to 
3500, and Pr is set to unity. 

Each run begins with a static, cylindrical magnetic flux tube embedded in a polytropic, fleld- 
free, reference state. The vertical domain of each simulation spans 5.147 pressure scale heights (or 
3.088 density scale heights). The tube initially has the form: 



B = VxVx^i-fVxJz 



(9) 



and 



Pqv = V X V X Wz + V y. Zz 



(10) 



(11) 



where 



Sx(r) 



(12) 



and 



Be{r) 



lrB,{r) . 



(13) 
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Here, BQ{r) denotes the azimuthal component of the field in the tube's cross-section, and Bx{r) 
refers to the axial component (which lies perpendicular to g = — g z along the Cartesian direction x 
). Both are given as functions of r, the radial distance to the central axis in the tube's cross-section. 
For each simulation, the initial size of the flux tube, a (deflned as the FWHM of Bx{r)), is set to 
O.lHr, and the magnetic field perturbations are scaled to the initial strength of the axial field at 
tube center (Bq). The distance over which a field line rotates once around the axis of the tube is 
given by 2TTa/q, where q is the non-dimensional twist parameter of equation (13). Note that this 
definition of q differs from that of Linton, Longcope, & Fisher (1996) by a factor of tube width, 
a. Both q and a are assumed constant, so that the initial rate of field line rotation per unit length 
along the tube (q/a for length scale a) remains fixed. To investigate how the tube's initial twist 
impacts the amount of fragmentation apparent during its rise, we consider three representative 
values of the twist parameter: q = 1/4, q = 3/16, and q = 1/8. 

Our goal is to model the dynamics of an emerging fi-loop. However, due to finite computational 
resources, we cannot evolve each flux tube self-consistently from an initial state of force balance (eg. 
Caligari, Moreno-Insertis, & Schtissler 1995; Fan & Fisher 1996; Fan 1999). We therefore introduce 
an ad-hoc, entropy perturbation at t = that causes the tube to rise and emerge in the shape of 
an r^-loop. We argue that the physics of the hydrodynamic interaction of the rising loop with its 
environment depends primarily on the geometry of the loop and its velocity fleld rather than how 
it arrived at its Jl-loop conflguration. The initial entropy perturbation is of the form 



Here, L refers to the initial length of the tube (which corresponds to the extent of the domain in the 
X direction), and S'o denotes the relative amplitude of the perturbation (taken to be unity in the 
dimensionless units of the code, where the unit of entropy is [s] = Cp[u]^/(r*[r])). With appropriate 
choices of L and length scale £, this initial condition has the effect of "pinning-down" the ends of 
the flux tube, while allowing the central portion to rise. For example, the run labeled SLl in Table 1 
has an initial acceleration due to buoyancy of 0.97 (in naturalized units of Bq / {8tt Hr Pr)) at the 
center of the tube, and a buoyancy contribution of —0.03 at each end. By varying the parameters 
in equation (14), we can investigate how the fragmentation process is affected by the geometry of 
a rising O-loop. 

Table 1 lists the q, L and £ parameter space that is explored in each of the 16 simulations, 
and assigns labels to each run. For convenience. Table 1 also shows the total number of field line 
rotations along the finite length of the tube, Q. Each horizontal flux tube is initially positioned 
near the bottom of the computational domain at zq = 0.1875 Zmax (-^max denotes the maximum 
vertical height of the domain). For the three-dimensional runs, the labeling convention consists 
of two letters ("L" or "S") followed by a number (0, 1, 2, or 3). The numbers denote the degree 
of magnetic field line rotation about the central axis of the tube (lower numbers imply a lesser 
amount of twist). The first letter of the label refers to the extent of the computational domain in 
the X direction (the "L" in this case refers to a "l"arge 512 zone domain, and the "S" stands for a 




(14) 
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"s" mall 256 zone domain) . The second letter of the label describes the eventual radius of curvature 
at the apex of the rising tube. If "L" is used, the run upon completion has a relatively "l"argc 
radius of curvature at the apex, and the fi-loop spans most of the computational domain in the x 
direction. Otherwise, "S" is used to denote a "s"horter O-loop; one which exhibits a smaller radius 
of curvature at its peak, and spans a smaller portion of the box length (this is accomplished by 
varying the parameter C in equation [14] in a fixed computational box). Labels with only one letter 
refer to the two-dimensional limiting cases. Note that it is the number of zones in the x direction 
that is changed between different cases, and not the grid resolution. 



3. Results 

In general, the magnetic field distribution in three dimensions can become quite complex. 
We find it useful to describe the field more intuitively in terms of the evolution and possible 
fragmentation of our initial magnetic flux tube. To accomplish this, we specify a new coordinate 
system based on the path of the flux tube through our Cartesian domain. We define the path 
of the tube in terms of the Cartesian variable x at any given time during a simulation. Setting 
up this new coordinate system and its basis vectors is a multi-step process. We first consider the 
magnetic field weighted moments of the position within vertical slices through relevant regions of 
the computational domain: 



z{x) = 



y{x) = -^-^ JJ y\B{x,y,z)\dydz, (15) 



where 



$(x) = jj \Bix,y,z)\dydz . (16) 
Then it is natural to define a path given by the vector 

ro{x) = XX + y{x)y + z{x)z (17) 
and to construct the Frenet tangent vector along this path. 

The path traced by tq is only a first approximation to the path of the flux tube. To establish 
the actual path of the tube, we calculate the magnetic field weighted position along a plane normal 
to £o passing through a given point along vq. This plane is defined by the equation £q- {r — fq) = 0, 
and its surface can be parameterized by r = u{v, w) x+v y +w z; where u{v, w) = {£y{y — v) +iz{z — 
w))/£x +x, and the 4's refer to the Cartesian components of Iq. With the area element along this 
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surface given by dS = \dr/dv x dr/dw\ dv dw = {k- £o) dv dw, the total unsigned magnetic flux 
across the plane can be expressed as $5 = Jg \'B{v,w) \ dS. We now define a new set of moments: 

V = — — // V \'B{v,w)\ {x • £o)~^ dv dw 



w = JJ w\B{v,w)\{x •£o)~'^dvdw. (19) 

Along with u = u{v,W), these points are used to define the path of a magnetic flux tube in a given 
region, 

r = ux + vy + WZ. (20) 

The frame of reference of the tube is then given by the Prenet tangent, normal, and binormal 
vectors along r: £ = dr/ds (where ds is the infinitesimal path length along r), h= K~^d£/ds , and 
b = £ X h , respectively. The tube's curvature at a given point along r is simply k = \d£/ds\. Of 
course, where the curvature of the flux tube is zero, the normal (and hence binormal) vectors are 
ill-defined. However, in the course of our analysis, we find that zero curvature occurs only at very 
localized n inflection points, or along horizontally oriented tubes (where vertical slices can be used 
to define a cross-sectional plane), and thus this limitation proves inconsequential. We note that 
our formalism is not terribly general, as it precludes the consideration of unusual configurations 
such as vertically oriented tubes, or tubes that are stacked on top of one another; however, none 
of these situations are encountered in our study. 

Depending on the initial degree of twist, portions of the rising tube will shed vortex pairs, 
and the magnetic flux will be redistributed in the tube cross-section such that much of the flux 
is located away from the tube's central axis. In some cases, the distribution of magnetic flux no 
longer resembles a single, cohesive tube; rather it has spread out and split apart into a configuration 
that can best be described as two separate flux tubes. At this point, we consider the tube to be 
"fragmented" . A quantitative measure of this fragmentation can be obtained by first calculating 
the second moments of the magnetic field strength along the normal and binormal directions of the 
tube: 

<{n-nf> = I {n-nf\B{n,b)\dS' 
^S' Js' 

<{h-hf> = ^ / (b-b)^\B{n,b)\dS' . (21) 
5>5' Js' 

Here, r = nh + bb spans the two-dimensional space of the plane normal to £ , dS' is the area element 
in that plane, and ^s' = fs' |B(^^, ^)| dS' is the corresponding unsigned flux (the £ dependence of 
the variables is implicitly understood). In the above equations, n = (l/^s') Js' ^ |B(?i, ^)| dS' and 
b = (1/<I>5') [c,, b |B(n, 6)1 dS' represent the first moments of the field distribution in this geometry. 
Due to the symmetry of our particular problem, the relative deformation of the tube along its 
cross-section can be represented by the ratio of the binormal to normal second moments: 

f= <b-W> _ (22) 
<(n-n)2> ^ ' 
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The total spread of the field distribution, cP' = {1/^s') Js' (('^ ~ '"')^ + {b — 6)^) |B(n, 6)| dS', can 
be used as a measure of how far the tube has dispersed in the cross-sectional plane. 

Prom an examination of numerous simulations, we empirically find that a magnetic flux tube 
can be considered to have split apart when its deformation, /, exceeds 1.5. We therefore use / 
as our measure of the degree of fragmentation of the tube, consistent with the results of Sehiissler 
(1979) and Longcope, Fisher, & Arendt (1996). If / > 1.5, then the path of each individual flux 
tube fragment is calculated using equations (19) and (20) over subdivided portions of the surface 
S. Since there is a fairly high degree of symmetry in our runs, the surface can be divided into two 
parts that lie on either side of a line defined by the Prenet normal. Thus, the regions of integration 
are easily determined. Figure 1 is a volume rendering of the magnetic field strength for the final 
timestep of run SLO, and shows the geometry of a tube which has fragmented as it rises toward 
the surface. The field strength distribution is visualized by casting parallel rays through the semi- 
transparent volume, and calculating the two-dimensional projection onto the viewing plane (IDL's 
"voxeLproj" routine). The red dotted line denotes the central axis of the tube through the volume. 
Where the tube has fragmented, yellow dotted lines trace the paths of each fragment. Note that 
the axis of each tube fragment does not coincide exactly with their centers of vorticity (where the 
field strength is locally concentrated); this offset depends on the amount of flux that resides in a 
thin "sheath" that extends from the forefront of the rising tube to the outer edge of each vortex 
pair. These features can be seen in a cross-section at the apex of this J7-loop in the bottom frame 
of Figure 2. In this case, the magnetic sheath is revealed as a relatively weak concentration of flux 
in the shape of a thin, semi-circular arch that lies above the primary concentrations of flux located 
in the trailing eddies. The yellow symbols denote the location of the central axis of each fragment, 
and the lengths of the blue dotted lines represent the second moments of the field distribution 
(given by equation [21]) at the centroid of the tube (note that at the apex of the loop, n and 
b approximately correspond to the Cartesian directions z and y respectively). The top panel of 
Figure 2 shows the apex cross-section of the flux tube shown in Figure 3 (the final timestep of run 
SS3). This is an example of a tube which, by our definition, has not yet fragmented as it approaches 
the photospheric boundary. 



3.1. How Fragmentation Depends on Tube Curvature and Twist 

To investigate how the degree of fragmentation at the apex of a rising, magnetic flux tube 
depends on its curvature at that point, we consider five separate runs: SSI, SLl, LSI, LLl, and 
LI. In each case, the initial twist parameter q is taken to be 1/8; that is, the simulations begin 
with identical, weakly-twisted horizontal flux tubes positioned at the same depth near the bottom 
of the computational domain. Effectively, the difference between each run is the length in the axial 
direction of the rising portion of the flux tubes. The changes in this length scale lead to substantive 
differences in the apex curvature of the i7-loop as it rises through the stratified atmosphere. Run 
LI (the two-dimensional limiting case) represents an infinitely long horizontal tube whose apex 
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curvature remains zero throughout its rise. Conversely, run SSI represents the shortest length 
scale we have considered, and thus refers to the run with the highest apex curvature as the tube 
approaches the photosphere. 

Figure 4 shows the degree of fragmentation along the axis of each tube for these five values of 
apex curvature at the time when each tube has risen approximately half-way through the vertical 
extent of the domain. In order of decreasing curvature, the simulations shown are runs SSI, SLl, 
LSI, LLl, and LI. The line / = 1.5 represents the (somewhat arbitrary) point at which we consider 
the tube to be fragmented. It is easy to see that as the level of apex curvature k, increases, the 
less the tube fragments for a fixed value of q. Figure 4 also shows that although the initial value 
of twist was slight, three of the four simulations have yet to show clear signs of fragmentation at 
the apex of the loop. This suggests that in three dimensions, the amount of twist necessary to 
prevent fragmentation is substantially reduced from the two-dimensional value — a point that we 
will return to later in this section. Note that toward the footpoints of each loop, the value of falls 
below 1 (the k, subscript is included to emphasize the dependence of the degree of fragmentation / 
on the tube's apex curvature k). This reflects the net elongation of the tube cross-section in the h 
direction as the il-loop expands. 

Figure 5 shows the time dependence of at the loop apex for the set of runs corresponding to 
those of Figure 4. A striking feature of this plot is the reduction by a factor of 1.5 in the degree 
of fragmentation near the photospheric boundary between the run with zero apex curvature, and 
the run with the maximum amount of curvature. The difference in is large enough that the 
magnetic flux emerges as two individual fragments when k = (run LI), yet is closer to being a 
single, cohesive tube when k > 1.7 x 10""^ (run SSI). Note that during the initial stages of the 
run, before the tube shows signs of fragmentation, f^^ falls slightly below 1. This is a result of the 
vertical elongation of the apex cross-section as the tube begins to rise, and flux is pulled into its 
wake. 

Wc find that if the initial value of the twist parameter q is greater than a "critical" value, gc, 
then the flux tube no longer emerges in a fragmented state, (/ < 1.5) regardless of the radius of 
curvature of the loop (see Figure 6). In our simulations, qc is empirically determined to be only 
slightly less than 3/16, the value chosen for the set of runs SS2, SL2, LS2, LL2, and L2 of Figure 7. 
For the two-dimensional limiting cases, the empirically determined critical value is consistent with 
the condition that the flux tube will split apart when the rise velocity of the tube exceeds the 
Alfven speed of the azimuthal field at its edge (note that the initial entropy perturbation in our 
simulations will affect the rise speed of the tube). However, if the initial twist of the flux tube is 
less than q^ we find that the curvature of the fi-loop plays an important role in determining the 
degree of fragmentation of the tube prior to its emergence through the photosphere — an effect not 
accounted for in previous two-dimensional studies of flux tube fragmentation (eg. Longcope, Fisher, 
&; Arendt 1996, Emonet & Moreno-Insertis 1998, and Fan, Zweibel, &: Lantz 1998). A comparison 
of Figure 7 with Figure 5 further suggests that as one chooses values of q progressively less than 
qc, the effect of loop curvature on the degree of apex fragmentation becomes more pronounced. 
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Run LI, the two-dimensional limiting case with a relatively small amount of initial field line 
twist, confirms the results of previous two-dimensional simulations (Longcope, Fisher, & Arendt 
1996; Emonet & Moreno- Inscrtis 1998; Fan, Zweibel, & Lantz 1998), which show that once frag- 
mentation has occurred, the two counter-rotating fragments repel one another via forces that result 
from infinitely long vortex lines. In two dimensions, the horizontal separation of the fragments 
can be understood in terms of the hydrodynamic force acting on an object with a net circulation 
r moving relative to a fluid with a velocity v: F = —p^vx T (see Fan, Zweibel, &; Lantz 1998). 
It is the component of this "lift" force per unit volume acting along the line between the centers 
of vorticity of the tube fragments that acts to push the tubes apart. In three dimensions, the 
fragments are finite in extent, and thus these forces are important only over a finite length of the 
loop. Figure 8 shows the trajectory of the two fragments (for the weakly-twisted case) at the tube 
apex for five different values of curvature. As the level of apex curvature of the tube increases, and 
the effective axial length scale of each vortex pair decreases, we see that the total volume integrated 
non-vertical component of the hydrodynamic lift — and thus the fragment separation — is reduced. 

Field line twisting due to non-uniform rotation of the fragments also acts to reduce the fragment 
separation. To illustrate the role of field line twist in the interaction of the vortex pairs, we consider a 
case where the azimuthal component of the magnetic field along the tube is initially zero. Figure 9 is 
a volume rendering of the magnitude of the current helicity density {\Hc\ = |J-B|, where J = V x B) 
for the final timestep of run SLO (where q ^0). Since the current helicity gives a measure of the twist 
of the magnetic field, Figure 9 shows that after the tube rises, spins differentially (see Figure 10), 
and splits apart, the magnetic field along each tube fragment becomes increasingly twisted. The 
magnetic forces imparted from the increased twist slow the rotational motion of the vortex pairs, 
and thus the tendency for the tubes to separate is further reduced. Figure 11 shows that as the 
apex curvature of the 0-loop increases, the net circulation near the apex of each fragment decreases. 
This lessens the repulsive force between fragments, and allows 0-loops with a high degree of apex 
curvature to behave more cohesively. This suppression of circulation in the tube fragments was 
predicted by Emonet & Moreno-Insertis (1998) and Moreno-Insertis (1997), who suggested that if 
the footpoint separation of an Jl-loop was small enough, the rotation of the vortex pairs would be 
suppressed. Since it is reasonable to assume a correlation between footpoint separation and apex 
curvature, we feel that the results of our simulations are generally consistent with this prediction. 



3.2. Field Morphology Prior to Emergence 

There arc a number of reasons why one must be cautious when comparing our results with flux 
emergence observations. First, the anelastic approximation becomes marginal as the magnetic flux 
tube approaches the photosphere, where densities decrease to the point that local sound speeds are 

comparable to Alfven speeds in the plasma. Additionally, in this first generation of models we do not 
include the effects of spherical geometry or the Coriolis force. Thus, predicted asymmetries in the 
geometry and field strength asymmetry of emerging active regions (see Moreno-Insertis, Schiissler, 
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& Caligari 1994 and Fan, Fisher &; DeLuca 1993) will not be reproduced. Further, we choose to use 
a simple, polytropic reference state rather than a more realistic convective background. As a result, 
the effects of convective turbulence on the tube arc neglected during its rise to the surface, likely 
over-simplifying the structure of the field. Nevertheless, we can make some general predictions 
regarding the morphology of emerging magnetic flux in active regions if the field exists in the form 
of a fragmented flux tube. 

Figures 12 and 13 make up a time-series of simulated vector magnetograms for the rising, 
fragmented flux tube of run SSI. These images were generated by taking a horizontal slice as 
close to the top of the computational domain as possible without having the topological evolution 
affected by interaction with the upper boundary. The grey-scale background represents the vertical 
component of the magnetic field; light regions correspond to outwardly directed magnetic field 
(toward the observer), and dark regions correspond to inwardly directed field (away from the 
observer). The arrows define the direction and relative strength of the transverse components of 
the field. Figures 14 and 15 show the corresponding velocity field at the same time intervals and 
locations. The first frame of Figure 12 shows the tip of the magnetic sheath as it begins to cut 
through the horizontal plane. At this point, there is little to no separation between regions of 
opposite polarity, and the corresponding flow field surrounding the magnetic region is relatively 
weak. This stage is short-lived, since the magnetic sheath is very thin. 

As the tube continues to rise, the horizontal cut passes through regions of magnetic field 
concentrated in the sheath on either side of the tube apex (which by now has emerged through 
the photosphere), as well as the field concentrated in the relatively horizontal trailing vortices of 
each fragment. As a result, an oval-like structure develops, as shown in frame 2 of Figure 12. The 
velocity field has a strong vertical component near the interfaces between the vortex tubes and the 
surrounding non-magnetized plasma, while exhibiting strong horizontal components at each end. 
The horizontal components reflect the flow of material along the sheath away from the apex (sec 
frame 2 of Figure 14). In frame 1 of Figure 13, the apex of each individual fragment has passed 
through the cutting plane, and we see only the portions of the tube fragments and sheath located 
just above the k inflection point of the O-loop. The two regions of opposite polarity now appear 
crescent-shaped as they continue to separate. After the fragmented portion of the loop has emerged 
through the photosphere, only the more concentrated, non-fragmented portions remain. Thus, in 
the flnal frame of Figure 13, the regions of strong magnetic fleld appear to have coalesced — a 
feature of emerging active regions that is commonly observed. 

Figures 16 and 17 show a similar time-series of vector magnetograms for runs LLl and SS3 
respectively. A comparison between these flgures and Figures 12 and 13 reveals the effect of initial 
twist and tube curvature on the overall characteristics of the emerging flux. A comparison of the 
surface separation of the individual fragments (the semi-minor axes of the ellipse-like structures 
shown in frame 2 of Figure 12 and frame 1 of Figure 16) can be used to determine the relative 
curvature between the two cases. In general, the lower the apex curvature of the fragmented fl- 
loop, the higher the "eccentricity" of the ellipse-like structure. Similarly, if the crescent-shape of the 



-13- 



opposite polarity magnetic structures is less distinct or absent in the simulations (compare frame 1 
of Figure 13 with frame 2 of Figure 17), this indicates that the initial level of twist of the magnetic 
flux tube was relatively high. Finally, asymmetries introduced by the higher level of twist present 
in the runs where g = 1/4 can result in the emergence of bipolar regions that arc tilted slightly 
with respect to the x direction. Note that a very slight tilt of approximately 4 degrees from the 
horizontal has developed between the center of the bipolar regions shown in frame 1 of Figure 17. 

Our description of a fragmented flux tube is in qualitative agreement with recent images of 
certain emerging active regions observed by Cauzzi, Canfield, k, Fisher (1996), Strous &; Zwaan 
(1999), and Mickey k, Labonte (1999). Vector magnetograms and white-light images of flux emer- 
gence taken with high temporal resolution by the Imaging Vector Magnetograph at the University 
of Hawaii's Mees Observatory (IVM) in relatively uncluttered portions of the solar disk (Mickey 
&: Labonte 1999) show evolving magnetic structures that are similar to the evolving structures 
described above. In particular, it is possible to infer the presence of oval-shaped magnetic features 
which evolve into crescent-shaped regions of opposite polarity that rapidly separate. In principle, 
direct measurements of the rate of separation and coalescence of the opposite polarity regions of 
the field may be used to obtain information about the structure and evolution of the sub-surface 
magnetic field. However, to perform this kind of detailed comparison with observation will require 
a more sophisticated set of models. 



4. Summary and Conclusions 

We have performed detailed numerical simulations of the rise of twisted magnetic flux tubes 
embedded in a non-magnetic, stratifled plasma using a code which solves the three-dimensional 
MHD equations in the anelastic approximation. The evolving magnetic field is described in terms 
of its volumetric flux distribution as a magnetic flux tube which may fragment into separate, distinct 
tubes during its rise toward the photospheric boundary. We find that the degree of fragmentation 
of the evolving magnetic fiux tube depends not only on the initial ratio of the azimuthal to axial 
components of the field along the tube (as was the case in two dimensions), but also on the three- 
dimensional geometry of the tube as it rises through the convection zone. The principal results of 
our analysis are the following: 

1. If the ratio of azimuthal to axial components of the field along a magnetic fiux tube exceeds a 
critical limit, then it will retain its cohesion and not fragment as it rises through the plasma. 
The critical limit occurs when the rise speed is approximately equal to the maximum Alfven 
speed of the azimuthal component of the field. This reinforces the conclusions of Emonet 
k, Moreno-Insertis (1998), Moreno-Insertis k Emonet (1996), Longcope, Fisher, k Arendt 
(1996), and Fan, Zweibel, k Lantz (1998). 

2. If the initial amount of fleld line twist is less than the critical value, and if the magnetic 
flux tube rises toward the photosphere as an O-loop, then the degree of apex fragmentation 
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depends on the curvature of the loop — the greater the apex curvature, the lesser the degree 
of fragmentation for a fixed amount of initial twist. 

3. In two dimensions, counter-rotating vortices are effectively infinite in extent and generate 
long-range flows which eventually prevent the continued vertical rise of tube fragments. This 
artificial geometric constraint is relaxed in three dimensions, and although the fragments 
continue to experience forces due to the interaction of the vortex pairs, those forces are only 
important along a short section of the tube. Thus, the fragments are able to rise to the 
surface. The forces due to vortex interaction depend upon the geometry of the O-loop — the 
greater the apex curvature, the lesser their magnitude. Thus, highly curved loops exhibit less 
fragmentation during their rise. 

4. Differential circulation between the apex and footpoint of an 0-loop leads to the introduction 
of new magnetic twist of opposite sign in each leg of a loop fragment. This twist reduces the 

circulation about the apex of each fragment, and further reduces the forces acting to separate 
the fragments of the flux tube. This result is consistent with the predictions of Emonet &; 
Moreno-Insertis (1998) and Moreno-Insertis (1997). 

5. Though these models do not admit to a direct, detailed comparison with observations, it is 
possible to infer certain general observational characteristics of emerging magnetic flux if the 
field configuration is that of a fragmented 0-loop rising through the solar surface. If the 
magnetic field erupts through a relatively quiet portion of the Sun, then we expect that one 
should observe concentrations of vertical flux which resemble expanding oval shapes. This 
type of feature should quickly evolve into longer-lived crescent-shaped regions of opposite 
polarity that steadily move away from one another. These extended regions should then 
coalesce once the fragmented portions of the r^-loop have emerged through the solar surface. 
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Fig. 1. — A volume rendering of for the last timestep of run SLO. The red dotted line denotes the 
central axis of the portion of the J7-loop which has not fragmented. The yellow dotted lines denote 
the paths of each distinct fragment in the region of the loop that is considered to be fragmented. 
The entire rectangular computational domain is shown. 
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Fig. 2. — The apex cross-sections of the il-loops shown in Figure 1 (bottom panel) and Figure 3 
(top panel). The symbols denote the positions of the tube axes, and the lengths of the blue dotted 
lines represent the second moments of the field distribution centered on the tube centroid (see 
equation [21]). These lines serve as a visual measure of the degree of fragmentation. The x and y 
scale is given in terms of computational zones. 



Fig. 3. — A volume rendering of for the last timestep of run SS3. The red dotted line denotes 
the central axis of the non-fragmented flux tube. 
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Fig. 4. — Fragmentation as a function of curvature for q = 1/8. The x scale is given by the number 
of computational zones from each tube's apex, and <z> is the average height of the central axis 
of the tubes (given in terms of the number of computational zones from the lower boundary). For 
clarity, the / values near the horizontal boundaries of each run (past the end of the Jl-loop) are 
not shown. 
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Fig. 5. — Apex fragmentation as a function of time and apex curvature for the set of runs where 
q = 1/8. 
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Fig. 6. — Apex fragmentation as a function of time and apex curvature for the set of runs where 
g = 1/4. 
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Fig. 7. — Apex fragmentation as a function of time and apex curvature for the set of runs where 
q = 3/16. 




Fig. 8. — Fragment trajectories at the apex of the rising J7-loop for the set of simulations where q = 
1/8. Note that at the loop apex, the normal n and binormal b directions correspond approximately 
to the Cartesian directions z and y respectively, ip represents the trajectory expected from pure 
expansion and no vortex interaction. The axis scale is in terms of computational zones — an entire 
vertical slice of the domain is shown. 



Fig. 9. — A volume rendering of the current helicity for the last timestep of run SLO. The dotted 
lines denote the paths of both the fragmented and un-fragmented portions of the loop. The entire 
rectangular computational domain is shown. 




Fig. 10. — The magnitude of the net circulation of an individual tube fragment (normalized to its 
value at the tube's apex, Tq) for the last timestep of run SLO. 
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Fig. 11. — The magnitude of the net circulation of a tube fragment (normalized to the value of 
run LI at t = 7, Fq) at the apex of the O-loops as a function of apex curvature for the runs with a 
relatively small amount of initial field line twist. 
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Fig. 12. — Vector magnetogram images for run SSI. The times given are in naturalized units 
{t = Hr/[v], see text), and the x and y scale are given in terms of the computational grid. The 
grey-scale background represents the vertical component of the magnetic field (positive values are 
indicated by the light regions), and the arrows represent the transverse components. 



Run SS 1 




Fig. 13. — Same as Figure 12, for later timesteps. 
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Fig. 14. — The velocity field for the same horizontal slice and times shown in Figure 12. The grey- 
scale background represents the vertical component of the velocity (light regions indicate positive 
values), and the arrows represent the transverse components of the velocity field. 
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Fig. 15. — Same as Figure 14 for the times corresponding to Figure 13. 
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Fig. 16. — Vector magnetogram images for run LLl. The times given are in naturalized units 
{t = Hr/[v], see text), and the x and y scale are given in terms of the computational grid. The 
grey-scale background represents the vertical component of the magnetic field (positive values are 
indicated by the light regions), and the arrows represent the transverse components. 
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Fig. 17. — Same as Figure 16 for run SS3. 



Table 1. Initial Flux Tube Parameters 
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SS2 


256x128x128 


4.362 


0.0852 


0.1875 


1.3017 


1.62260 


SS3 


256x128x128 


4.362 


0.0852 


0.2500 


1.7356 


1.62906 


SLO 


256x128x128 


4.362 


0.1704 


0.0000 


0.0000 


1.05930 


SLl 


256x128x128 


4.362 


0.1704 


0.1250 


0.8678 


1.04928 


SL2 


256x128x128 


4.362 


0.1704 


0.1875 


1.3017 


0.99512 


SL3 


256x128x128 
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0.1704 
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1.7356 


1.02579 


LSI 


512x128x128 


8.724 


0.1704 


0.1250 


1.7356 


0.53913 


LS2 


512x128x128 


8.724 


0.1704 


0.1875 


2.6034 


0.51325 


LS3 


512x128x128 


8.724 


0.1704 


0.2500 


3.4712 


0.46269 


LLl 
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0.2500 


oo 


0.00000 



^in units of Hj. 

■^apex curvature for last timestep of each run (xlO""^) 
^two-dimensional limiting case 



